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Abstract 

Extensive Monte-Carlo simulations were performed to evaluate the excess 
number of clusters and the crossing probability function for three-dimensional 
percolation on the simple cubic (s.c), face-centered cubic (f.c.o), and body- 
centered cubic (b.c.c.) lattices. Systems L x L x L with L » L were 
studied for both bond (s.c, f.c.c, b.c.c.) and site (f.c.o) percolation. The 
excess number of clusters b per unit length was confirmed to be a universal 
quantity with a value b ~ 0.412. Likewise, the critical crossing probability in 
the L' direction, with periodic boundary conditions in the L x L plane, was 
found to follow a universal exponential decay as a function of r = L /L for 
large r. Simulations were also carried out to find new precise values of the 
critical thresholds for site percolation on the f.c.c. and b.c.c. lattices, yielding 
p c (f.c.c.) = 0.199 236 5 ±0.000 0010, p c (b.c.c.) = 0.245 961 5 ±0.000 001 0. We 

also report the value p c (s.c.) = 0.311 608 ± 0.000 000 4 for site percolation. 
PACS numbers(s): 64.60Ak, 05.70.Jk 
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I. INTRODUCTION 



The standard percolation model 0] involves the random occupation of sites or bonds of a 
regular lattice. At a critical occupation probability p c , the mean size of clusters of occupied 
sites becomes infinite, while the number of clusters n(p) per site or per unit volume remains 
finite with n c = n(p c ). 

The value of n c depends on the microscopic characteristics of each system, and because 
of this it is a non-universal quantity For two-dimensional (2d) systems, precise numerical 
values of n c for bond and site percolation on the square and triangular lattices were found by 
Ziff, Finch, and Adamchik 0, whose results for bond percolation confirmed the theoretical 
predictions of Temperley and Lieb |§ and Baxter, Temperley, and Ashley @] . In 3d, there 
are no theoretical predictions for n c , and its values for different systems apparently have not 
been reported in literature. 

In Ref. [[|], it was also found that the excess number of clusters b = lim^oo LL' (n(L, L') — 
n c ), with r = L /L = fixed, where n(L, L') is the number of clusters per unit area in a finite 
system of size LxL' with periodic boundary conditions, is a universal quantity that depends 
only upon aspect ratio r. (Note that in 0], the authors defined n as clusters per site rather 
than per unit area, but the result for b is the same.) This universality is consistent with the 
arguments of Privman and Fisher []5| , and has also been discussed by Aharony and Stauffer 
H and by Miiller 0] for the Ising model. Kleban and Ziff || introduced an excess number 
per unit length b = lim^oo b(r)/r = lini^oo L 2 (n(L, L') — n c ) for long cylindrical systems 
L' >> L, and derived exact results for both b(r) and b in 2d systems. Again, however, no 
theoretical predictions for b in 3d exist. 

In this paper, we determine n c and b for various 3d rectangular solid systems of dimen- 
sions L x L x L' with L' >> L. We consider bond percolation on the simple cubic (s.c), 
body-centered cubic (b.c.c), and face-centered cubic (f.c.c.) lattices, and site percolation 
on the f.c.c. lattice. 

A prerequisite to finding the value of n c for each of these systems is knowing the critical 
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occupational probability p c to high accuracy. Previously, accurate values were found for bond 
percolation on all three lattices and site percolation on the s.c. lattice only, as summarized 
in Table |. To round out these values, we carried out simulations to determine p c for site 
percolation on the b.c.c. and f.c.c. lattices to high accuracy — although we used only the 
latter in the study of the excess cluster number, since the universality was clearly confirmed 
with the four systems that we studied. In another work U, we have studied site percolation 
on the s. c. lattice, and report this result in Table | also. 

The simulations for finding n c were also used to study the critical crossing probability for 
the three-dimensional systems. The crossing probability function 7r(r) of a system of shape 
T gives the probability that at least one cluster connects two disjoint pieces of the boundary 
dT, and has been of much interest lately following the realization that it is a fundamental, 
universal property of percolation, independent of the underlying lattice type, and subject 
to conformal invariance [pTDj p~3|] . In 2d, Cardy [|ll]] derived an explicit expression for the 



vertical crossing probability n v of rectangular systems L x L , with open boundaries in the 
horizontal direction, and Watts [TH| derived an expression for the probability of vertical but 



not horizontal crossing for this system. The ir v for 2d systems with periodic (and other) 
boundary conditions was studied by Hovi and Aharony |T5|]. A number of systems were 
also studied by various groups including Hu et al. [[U|, Hsu et al. |l7j, Gropengiesser |l8 



and Vicsek and Kertesz [fil| . In 3d, work has been restricted to simple cubical boundaries 



LxLxL, with crossing studied between two opposite planes and various boundary conditions 



on the sides [16,20,21 



Here we find ir v for the L x L x L' systems for all L' by measuring the distribution 
of the maximum height of clusters connected to the base of the rectangular system. (A 
similar method was used in |22] for 2d systems.) We consider crossing in the L' direction 
for systems with periodic boundary conditions in the L x L plane, and show that tt v is a 
universal function of r = L' / L for large L. 

In the following three sections we report on the determination of the new values of p c , 
the determination of n r and 6, and the determination of 7r„(r). The results are summarized 



and discussed further in the conclusions section. 



II. PERCOLATION THRESHOLDS 

Precise values for the thresholds for bond percolation on all three lattices, and for site 
percolation on the s.c. lattice, have been found elsewhere. Here we also determine accurate 
values for site percolation on the f.c.c. and b.c.c. lattices. A summary of our results and 
other recent results is given in Table |. 

The procedure we used to find p c was similar to that we used for bond percolation in 
23] . We grew individual clusters by a Leath-type algorithm and identified the critical point 



using an epidemic scaling analysis. A virtual lattice of 2048 3 sites was simulated, using 



the block-data method first described in [24J]. There were only two minor changes made 



to the simulation of [^3| so that it could be used to study site percolation. First, as the 
clusters were grown, the sites were either occupied with a probability, p, or left vacant with 
a probability, 1 — p. If a site was determined to be vacant, then (unlike in bond percolation) 
it was never revisited as a potential growing site. The other difference is the cut-off for the 
growth of these clusters was set to 2 19 (524,288) wetted sites, as opposed to 2 20 (1,048,576) 



and 2 21 (2,097,152) in ref. |23|. 

The simulation yielded the fraction of clusters P(s,p) that grew to a size greater than 
or equal to s sites. When p is near p c , one expects P(s,p) to behave as 

P{s, p) ~ As 2 - T f((p - p c )s°) « As 2 ^[l + C{p - p c )s° + ...] (1) 



where r and a are universal exponents [25]. Here we assumed the values r = 2.189 and 



a = 0.445, consistent with other 3d work [f23|, p6| , p7|j . As in |23|, plots of s T ~ 2 P(s,p) vs. 
s a for site percolation of the b.c.c. and f.c.c. lattices were used to find the value of the 
percolation threshold, which corresponds to horizontal behavior for large L on such a plot. 
In all, we generated 1.5 x 10 7 clusters for the f.c.c. lattice and 2.2 x 10 7 for the b.c.c. lattice 
for a range of values of p requiring several weeks of workstation computer time. The results 
are plotted in Fig. [l| and imply the following values for the critical thresholds: 



p c (b.c.c.) = 0.245 961 5 ± 0.000 001 

p c (f.c.c.) = 0.199 236 5 ± 0.000 001 (2) 

These results were consistent with (and more than 1000 times more precise than) previous 
work, as shown in Table |. 

III. VALUES OF N c AND THE FINITE-SIZE CORRECTION B 

Using the values of the critical thresholds given in Table |, we carried out simulations 
to measure the number of clusters for bond percolation on each of the three-dimensional 
lattices and site percolation on the f.c.c. lattice. (We did not consider site percolation on 
the s.c. and b.c.c. lattice in this calculation.) Clusters were grown successively from every 



unvisited site by a growth algorithm |28j on a three-dimensional square bar, L x L x L 
with L' » L. Periodic boundary conditions were assumed in each horizontal plane. (Here, 
vertical is taken to be the L' direction). The first cluster was started in the upper left- 
hand corner of the first plane (z = 0) at the point (0,0,0). From this corner, a cluster was 
grown to the nearest neighboring sites as defined for each system by the unit vectors in |23 



occupying the connecting bonds or neighboring sites with a probability, p c , and leaving them 
unoccupied with a probability, 1 — p c . After the first cluster was grown, a new cluster was 
seeded from the first unoccupied site in the left-most column, and grown until it died. After 
all sites of the first plane were tested, the growing plane was moved to z = 1, and so on. 
Because the previous planes were completely occupied, their data could be discarded and 
the memory recycled. Furthermore, the clusters never extended up to a plane of distance 
z = 32L from the growing plane. As a consequence, a system of size L x L x 32L could be 
used to effectively simulate aLxixoo system by wrapping around in the third direction. 

We ran simulations to L' = 2 x 10 9 , with L — 4, 5, 6, 7, 8, 10, and 12 for the s.c. 
lattice, L = 4, 6, 8, 10 and 12 for the f.c.c. lattice (both site and bond), and L — 4, 6, 8, 
10, 12, and 16 for the b.c.c. lattice. In total, we grew about 1.06 x 10 12 clusters for the 
b.c.c. lattice, 1.88 x 10 12 clusters for the s.c. lattice, 1.44 x 10 12 clusters for bond percolation 



on the f.c.c. lattice, and 3.32 x 10 11 for site percolation on the f.c.c. lattice, which required 
several additional months of computer time. 

In Figure 0, we display a representative 4x4 plane of each of the three lattices, showing 
how the lattices were oriented in our simulations and how the unit dimension was defined. 
The darkened circles represent active sites in the current plane, and the empty circles rep- 
resent active sites that are in the neighboring planes. The solid lines are bonds which lie 
within the current plane and the dashed lines represent bonds which connect the displayed 
plane to the neighboring planes. For modeling the s.c. lattice, the plane shown in Figure 
H](a) is repeated for the whole length of the cylinder, while for the other two lattices, the 
plane shown in the figure is repeated on every other plane. In the case of the s.c. lattice, 
all of the available sites in the plane are considered active, for the f.c.c. lattice, only half of 
the underlying cubic-lattice sites are active, and for the b.c.c. lattice, only a quarter of the 
cubic-lattice sites are active. Note that the unit dimension that we define for the b.c.c. and 
f.c.c. lattices is neither the unit cell dimension nor the nearest-neighbor distance, but half 
of the unit cell dimension. 

Now, for a finite system of volume V with periodic boundary conditions, analogous to 
what was found in [g] for 2d, we expect 

b c 

n = ric+ V + V 2+ '" ™ 

where b, representing the excess number of clusters in this finite system, is universal, a 
function of the shape only. Here we studied L x L x L' systems, where L' >> L, with the 
volume given by V = L 2 L' . For systems of this shape, we expect the excess number of 
clusters per unit length b/(L / L) to be a constant 6, i.e., 

b ~ bL'/L (4) 

for L' >> L. Likewise, we write c ~ c(L' /L) 2 . Then it follows from ([|) that 

n = n c + — + — + .. . (5) 
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Both b and b are functions of the system shape only and are universal quantities, but c and 
c vary from system to system and are not universal. Eq. fl5|) implies that n c can be found 
from a plot of n vs. 1/L 3 , as shown in Figure |3| for our data from the b.c.c. lattice. The 
values of n c , which are shown in units of number of clusters per unit volume as defined in 
Fig. for the various lattices, are given in Table [H]. These values can be converted to units 
of number of clusters per site by taking into account that the s.c, f.c.c, and b.c.c. lattices 
have 1, 1/2, and 1/4 sites per unit volume, respectively. 
Equation (fj|) can be rearranged as 

(n - n c )L 3 = b + + . . . (6) 

Therefore, once n c is determined, b and c can be found from a plot of (n — n c )L 3 vs. 1/L 3 . 
Figure f| shows this plot for the systems that we studied. The resulting values of b and c for 
each of the systems are listed in Table |J. A universal value of b = 0.412 ± 0.002 is obtained 
from these results. 



IV. CRITICAL CROSSING PROBABILITY vr y 

Our simulations for n c could also be used to obtain ir v by comparing the distance from 
the growth plane to the maximum height plane. If this distance is greater than or equal to 
some fixed value L' , then crossing will occur in an L x L x L system (with periodic boundary 
conditions in each L x L plane). In other words, we could determine ir v (L, L,L) for all L' 
by keeping track of the distribution of distances between the growth plane and maximum 
height planes in our continuous simulations. 

In 2d, the probability of crossing a system of aspect ratio r = height/width in the vertical 



direction, with periodic boundary conditions in the horizontal direction, is given by ||13| , [29 

7r„(r) ~ e - 2 " r(2 - D) = e~^ r (7) 

for large r, where D = 91/48 is the 2d fractal dimension. Eq. (|7|) follows from a conformal 
transformation from an annulus to a rectangle, using that the probability a cluster extends 



beyond a radial distance R scales as R D ~ d . We have separately verified that Eq. ([?]) holds 
accurately for all r somewhat greater than 1. 

For 3d systems, while it still is true that the radial probability scales as R D ~ d , we 
cannot connect it to n v of the L x L x L' system, because we cannot make a conformal 
transformation between the concentric spheres and a rectangular solid. However, we still 
expect an exponential dependence upon r = L' /L, because that term represents the smallest 
eigenvalue of the transfer matrix. We thus hypothesize 

ir v ~ Ke~ mr . (8) 

for large r. To check this, we plot hnr v vs. r in Fig. ^j, which contains the results from 
all four systems studied, for L — 8, 10, and 12. To get the best data collapse, we defined 
r = (L' +£)/L, which allows for a lattice finite-size effect or boundary extrapolation length in 
the L' direction, in which the effective location of the free boundary is not uniquely defined 
f2~2fl . (Such an ambiguity in size does not occur in the L directions, because of the periodic 
boundary conditions.) In fact, the data for all three bond percolation systems collapsed 
nicely with i = —1.3, while the data for site percolation on the f.c.c. lattice required a 
constant of I = 1.36 to fall on the same curve. Figure § shows the effect of I by comparing 
an enlarged portion of our data from the s.c. (bond) lattice when i = and I = —1.3. The 
corresponding values of m and hiK are —1.37 ± 0.01 and 0.75 ± 0.05, respectively. 

V. DISCUSSION OF RESULTS 

Our values for the critical thresholds of site percolation on the f.c.c. and b.c.c. lattices 
are listed in Table |. Along with the other results which are summarized in that table, the 
thresholds of all three 3d systems, for both site and bond percolation, are now known to a 
very high accuracy. 

Table lists n c , b and c for the four systems studied. Our simulations confirm that b is 
universal in 3d as it is in 2d 0], with a value b ~ 0.412. In 2d, the corresponding value is 
6 = 5^3/24 = 0.360 844... @. 
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The average density of clusters per site, n c , varies from system to system, as expected. 
The values for n c in Table [TI| show that the simple cubic is the most dense system, according 
to the convention we used to define the unit volume of the system. 

Our simulations have also shown that n v is universal as shown in Figure ^j, and 
possesses an exponential decay (Q) with m = 1.37 ± 0.01, compared with a value of 
57r/24 = 0.654498... in 2d. For a cubical system (L x L x L or r = 1), Eq. (||) im- 
plies a value of tt v = 0.54 ± 0.04, while a direct analysis of our data at that point yields 
the more precise value tt v = 0.573 ± 0.005. The latter value is somewhat higher than the 
result 0.513 ± 0.005 recently reported by Acharyya and Stauffer |ZT| for a system with heli- 
cal boundary conditions in the plane, which are similar to periodic boundary conditions but 
with the rows shifted by one. We believe that in the limit of large L these two boundary 
conditions should be equivalent, although this belief is not supported by the discrepancy in 
the values seen above. 

Many additional questions are raised for 3d systems. What is b(r',r") where r' = L' /L 
and r" = L" / L for an L x L' x L" system (with periodic boundary conditions in all direc- 
tions)? What is the effect of helicity or a twist of the order L in the periodic boundary 
conditions? Is b related to the number of "percolating" clusters per unit length (however 
precisely that may be defined)? Finally, can one devise a system that conformally transforms 
to concentric spheres, so that the crossing probability across that system will be given by a 
formula analogous to (^)? 
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FIGURES 

FIG. 1. Plot of s T ~ 2 P(s,p) vs. s a for (a) f.c.c. and (b) b.c.c. lattices using r = 2.189 and 
a = 0.445. Each curve represents a different value of p, which are (from top to bottom) (a) 
0.199 2375, 0.199 236 5, and 0.199 235 5, and (b) 0.245 962 5, 0.245 9615, and 0.245 960 5. 

FIG. 2. Representative 4x4 planes for the (a) s.c, (b) f.c.c, and (c) b.c.c. lattices. The 
darkened circles represent active sites in the plane and empty circles represent active sites in the 
neighboring planes. The solid lines represent bonds in the plane and dashed lines represent bonds 
which go to the neighboring planes. 

FIG. 3. Plot of n vs. 1/L 3 for bond percolation on the b.c.c. lattice. The intercept of this 
plot yields n c and the slope yields b according to Eq. (g). 

FIG. 4. Plot of (n — n c )L 3 vs. 1/L 3 for the b.c.c. (bond), f.c.c. (bond), f.c.c. (site), and s.c. 
(bond) (from top to bottom) systems at p c . In these plots, the intercept represents the value of b 
and the slope is the second correction term c. The values of b and c are listed in Table [J 

FIG. 5. Plot of lmr v vs. r = (L +£)/L for the b.c.c. (bond) (dashed lines), f.c.c. (bond) (dotted 
lines), f.c.c. (site) (also dotted lines), and s.c. (bond) (solid lines) lattices of size L x L x L with 
L = 8, 10, and 12 at p c - In these plots, the intercept represents the value of InK and the slope is 
m. The values of \n.K and m are 0.75 ± 0.05 and —1.37 ± 0.01. 

FIG. 6. Plot of hi 7^ vs. r = (l! + t)/L for an enlarged portion of the data from the s.c. (bond) 
lattice of size L x L x L' with L = 8 (square), 10 (circle), and 12 (triangle) at p c . The upper 
three curves show the data plotted with I = 0, and the bottom curve shows the data collapse when 
t = —1.3 is used. 
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TABLES 

TABLE I. Values of p c for bond and site percolation on the b.c.c., f.c.c. and s.c. lattices from 
present (*) and other recent work. The numbers in parenthesis represent the errors in the last 
digit (s). 



system 



b.c.c. (bond) 



b.c.c. (site) 



Pc 



0.180 3 
0.180 2(2) 
0.180 287 5(10) 
0.246 
0.245 8(2) 
0.245 9615(10) 



Ref. 



m 



Value used here 



0.180 2875 



f.c.c. (bond) 



f.c.c. (site) 



0.119 
0.1200(2) 
0.120163 5(10) 
0.198 
0.199 4(2) 
0.199 236 5(10) 



23 



m 



0.120163 5 



0.199 236 5 



s.c. (bond) 



s.c. (site) 



0.248 8 
0.248 7(2) 
0.248 8(2) 
0.248 75(13) 
0.248 814(3) 
0.248 812(2) 
0.248 812 6(5) 
0.3116 
0.3114(4) 
0.311605(10) 



m 
m 

m 
m 

1231 



30 
26 



0.248 812 6 
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0.311604(6) | 
0.311605(5) | 
0.311600(5) | 
0.3116081(13) | 
0.311608 0(4) | 


a 

21 1 
34| 
27| 

1 


TABLE II. 


Values of n c (clusters per unit volume), b, and c for the systems 


studied. 


system 


n c 


b 


c 


s.c. (bond) 


0.272 9310(5) 


0.414(3) 


6.0(7) 


f.c.c. (bond) 
f.c.c. (site) 


0.153 844 0(5) 
0.013 265 5(5) 


0.414(3) 
0.409(3) 


-1.4(3) 
-1.8(3) 


b.c.c. (bond) 


0.074 586 0(5) 


0.409(3) 


-5.5(7) 
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Figure 1(b) 
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